The present invention relates to seismic migration and particularly its use in imaging and, more particularly, but not exclusively to seismic migration using a Krylov space expansion of the square root operator to approximate a wave equation to allow rapid and accurate modeling of wave propagation between layers.
Geological surveying has widespread applications, for example to locate oil and gas reserves whether on land or offshore. In the past, surveying was carried out by looking at surface geological formations and using the surveyor""s knowledge and experience to determine locations of underground or subsurface structures likely to contain reserves. A trial bore was then made at the determined location and tested for the presence of hydrocarbons.
The above process owed much to trial and error, and trial bores are expensive. There is therefore a need to reduce the number of trial bores needed for successful discovery. Thus, more recently methodology has been developed to assist the surveyor by allowing him to learn about the sub-surface structure during the initial survey and prior to making a trial bore. Knowing about the substructure allows for a more educated location of the trial bore and thus improves the efficiency of the surveying process.
One method of finding substructure formations uses satellite imaging. Another method uses seismic imaging. In seismic imaging, a technique similar in principal to radar and sonar is used, in which the time taken by sound waves to travel through the sub-surface structure and return to the surface is measured and used to infer the structure. Seismic data is obtained by sending an energy pulse (basically a sound wave) into the earth and then listening for it to be reflected off rock layers and return to the surface. The time it takes for the pulse to return indicates how far it has traveled. The direction and timing of the received wave indicate where it bounced off a rock layer and allows an image to be constructed indicating the position of those layers.
Seismic imaging can be used in conjunction with software analysis to produce three dimensional images of subsurface rock formations. Seismic imaging thus reduces the trial and error dimension of surveying and provides for greater efficiency.
In more detail, sound waves are produced by a seismic source. The seismic source may for example be a small underground explosion referred to hereinafter as a shot. Dynamite is common on land, and air guns which produce large bubbles are common in water. Another way of producing sound waves is known as Vibroseis, in which a heavy vehicle is shaken in such a way as to produce a set of vibrations. The sound waves subsequently propagate into the earth and partially reflect from interfaces, across which the subsurface velocity or the density varies discontinuously. The reflected waves are recorded by recording devices, known as geophones (or hydrophones in marine surveys), which are usually placed in the vicinity of the surface, or within well bores. The time history of each device, which records the amplitude of the waves detected following a given shot, is stored as a seismic trace. The objective is to deduce the sub-surface structure from the recorded data.
Important in constructing an image from shot data are the relative locations of shot sources and geophones as well as information of the actual energy source used. The relative locations of shot source and geophones dictates the region of the subsurface that is imaged and the resolution level of that image.
Seismic migration or imaging is the process by which the seismic data is mapped to form an image of the subsurface. In general, this mapping requires knowledge of the subsurface velocity. The subsurface velocity is highly variable, depending on the type of material in the subsurface structure. The fact of such velocity variability in the first few kilometers beneath the Earth""s surface makes the task of seismic imaging more difficult than other types of acoustical imaging, such as for example ultra-sound and sonar which are carried out through a broadly unified medium. A description of the seismic method and of seismic imaging can be found in standard texts such as M. Dobrin and C. Savit: Introduction to Geophysical Prospecting. McGraw-Hill, New York. 1988, O. Yilmaz. Seismic Data Processing. Society of Exploration Geophysics. 1987, the contents of which are hereby incorporated by reference.
Depth migration refers to a type of imaging which maps the time recorded input seismic traces, that is the time history recorded at each device referred to above, into a subsurface spatial image, taking into account the different velocities at different depths in the structure. Migration is based on solution of a governing wave equation. Most often the acoustic wave equation is used because of its simplicity, and because it produces the correct arrival times of primary (P) waves. However, when more accurate results are required, the elastic equations, or other more elaborate equations may be used. The specific embodiments of the imaging technique disclosed herein use the acoustic wave equation. However the same technique can be applied to imaging based on other types of wave equations, such as the elastic wave equations.
Depth migration methods can be roughly divided into two categories; methods based on direct solution of the wave equation, often termed wave equation migration, and methods based on geometrical optics based approximate solutions to the wave equations. The latter have been termed in the exploration industry as Kirchhoff migration. Kirchhoff migration has been considered faster and more flexible in accepting irregular source-receiver geometry. The downward propagation of the surface data as modeled in wave equation migration is not based on geometrical optics approximations and is thus more realistic although harder to compute. Consequently wave equation migration type approaches are considered to have the potential of producing more accurate results. The embodiments of the present invention described below apply to wave equation migration.
The Basis of Wave Equation Imaging
Considering current art wave equation imaging in greater detail, the surface recorded seismic data can be grouped in different ways, for example, according to shots. In grouping according to shot, all seismic traces produced by a given shot are aggregated together into what is known as a common shot gather. As an alternative, seismic traces may be grouped according to receivers. That is all traces recorded by a surface receiver are aggregated together into a common receiver gather. As a third possibility, grouping according to offsets, all traces for which the shot-receiver separation falls within a specified range are aggregated together into a common offset gather. Each of the above types of data grouping, and other, similar possibilities, has its own imaging method, but the general imaging principles behind all methods are similar. For the sake of clarity and conciseness, the following discussion is limited to acoustic common shot imaging of surface recoded reflection data. The skilled person will appreciate how to apply the principles to the other methods.
Let p(x,y,z=0, t) represent the time history of receivers for a given shot, where (x, y) are the horizontal coordinates of each receiver position. The plane (z=0) defines the surface of the earth, which in the present discussion is assumed to be flat. If needed the scheme can be modified to account for topography. The seismic migration, or imaging, consists of two steps:
1. Extrapolation of the surface data p(x,y,z=0,t) in depth to form p(x,y,z,t).
2. Application of an imaging condition to create the subsurface image Pmig(x,y,z).
The subsurface image can be calculated according to pmig(x,y,z)=p(x,y,z,t=td), where td is the time of arrival of a direct wave from the shot location to the subsurface location, see M. Reshef and D. Kosloff. Migration of common shot gathers. Geophysics 51, 321-331. 1986, the contents of which are hereby incorporated by reference. An alternative imaging condition which is often used is obtained by cross correlating, Pmig(x,y,z)=p(x,y,z,t)*Pmodel(x,y,z,t)|t=0, where Pmodel(x,y,z,t) denotes the calculated wave amplitude which results from a point source located at the shot position. There are other variants of the imaging condition, all of which use the extrapolated wave amplitude.
The main differences between the imaging techniques available in the seismic industry lie in the method of extrapolating the surface data. The extrapolations are based on solution of the wave equation. For an acoustic medium with constant density, the wave equation reads,                                           1                          c              2                                ⁢                                                    ∂                2                            ⁢              p                                      ∂                              t                2                                                    =                                                            ∂                2                            ⁢              p                                      ∂                              x                2                                              +                                                    ∂                2                            ⁢              p                                      ∂                              y                2                                              +                                                                      ∂                  2                                ⁢                p                                            ∂                                  z                  2                                                      .                                              (        2.1        )            
t denotes time, and c(x,y,z) is the seismic velocity. An assumption of constant density is often used since it produces the same arrival times as the variable density equation and it is simpler to solve. An alternative expression is obtained after a temporal Fourier transform of this equation,                                                         ∂              2                        ⁢                          p              ~                                            ∂                          z              2                                      =                                            -                                                ω                  2                                                  c                  2                                                      ⁢                          p              ~                                -                                                    ∂                2                            ⁢                              p                ~                                                    ∂                              x                2                                              -                                                                      ∂                  2                                ⁢                                  p                  ~                                                            ∂                                  y                  2                                                      .                                              (        2.2        )            
{tilde over (p)}(x,y,z,xcfx89) denotes the Fourier transform of p(x,y,z,t), and xcfx89 denotes the temporal frequency. Equation (2.1) is most suitable for extrapolation in time. The resulting type of migration is called reverse time migration, see G. McMechan. Migration by Extrapolation of Time Dependent Boundary Values. Geophysical Prospecting, 31, 413-420. 1983; and E. Baysal, D. Kosloff, and J. Sherwood. Reverse Time Migration. Geophysics, 48, 1414-1524, 1983. The contents of both of these citations are hereby incorporated by reference. With the reverse time migration approach, the surface recorded data is introduced as a surface boundary condition for solving (2.1), where the values are introduced in a reverse order in time. The reverse time type of migration has the potential for high accuracy, however it is considered to be slow. A more efficient type of migration is obtained by a downward extrapolation in depth based on (2.2). Accordingly, each frequency slice {tilde over (p)}(x,y,z=0,xcfx89) is used as an initial condition for calculating {tilde over (p)}(x,y,z,xcfx89) at all depth levels. These methods are often termed fxy migration methods, where ƒ stands for frequency. The problem with downward extrapolation is that in the case of variable velocity, in fact the usual case, it is not possible to solve (2.2) by a direct integration in z. This can be understood by examining the solution of (2.2) in the wave number domain when the velocity c is constant,                                                         p                              ~                                  ~                  ~                                                      ⁡                          (                                                k                  x                                ,                                  k                  y                                ,                z                ,                ω                            )                                =                                                    A                +                            ⁢                              ⅇ                                  ⅈ                  ⁢                                      xe2x80x83                                    ⁢                                      k                    2                                                                        +                                          A                -                            ⁢                              ⅇ                                                      -                    ⅈ                                    ⁢                                      xe2x80x83                                    ⁢                                      k                    2                                                                                      ,                            (        2.3        )            
where       p          ~              ~        ~              ⁡      (                  k        x            ,              k        y            ,      z      ,      ω        )  
is the two dimensional spatial Fourier transform of {tilde over (p)}(x,y,z,xcfx89), and             k      z        =                                                      ω              2                                      c              2                                -                      k            x            2                    -                      k            y            2                              ·              k        x              ,      k    y  
are respectively the x and y wave numbers. When                     k        x        2            +              k        y        2               less than                   ω        2                    c        2              ,      k    2  
is real, the solution represents propagating waves. In this case A+ and Axe2x88x92 represent the amplitudes of down-going and up-going waves respectively. Conversely when                     k        x        2            +              k        y        2               greater than                   ω        2                    c        2              ,      k    z  
is imaginary and the solution represents exponentially decaying or increasing solutions, depending on the sign in the exponential. The exponentially increasing solution components cause numerical instability and therefore need to be eliminated. While this elimination can be easily done in (2.3), in the case of general velocity variation c(x,y,z), these components cannot be isolated in a simple manner. Existing methods for fxy migration attack this difficulty in a number of ways.
One approach for downward continuation uses one-way wave equations instead of the acoustic wave equation, and reference is made to J. Claerbout and S. Doherty. Downward continuation of move-out corrected seismograms: Geophysics 37, 112-139. 1972 the contents of which are hereby incorporated by reference. The solution to these equations contains only up (or down) going propagating waves without exponentially growing components. In regard to the propagating wave components, the equations have solutions which are close in value to the solutions to the exact wave equation in cases up to a certain dip angle. The limiting dip angle depends on the complexity of the particular one-way wave equation being used, where as a rule, a better and steeper dip response requires a more complex equation and a larger numerical effort to solve it. The main drawbacks of this approach are numerical dispersion and a dip angle limitation. Consequently, migration approaches based on one-way wave equations may fail to image structures containing steep dips.
A second class of downward propagation methods is based on a design of explicit operators for a suite of constant velocities, and reference is made to D. Hale. Stable explicit depth migration of seismic wavefields. Geophysics 56 1770-1777, 1991; and R. Soubaras. Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis. Geophysics 61, 1386-1393. 1996, the contents of both of these documents being hereby incorporated by reference. The downward continuation from one depth level to the next becomes an application of a spatially variant filter, where the filter coefficients at each image point are selected according to the value of the velocity thereat. The explicit operators approach also suffers from dip limitations since the filter coefficients are designed to match the solution (2.3) up to a specified dip value. It is noted that for a given frequency xcfx89 and constant velocity c, the dip angle xcex8 is related to the horizontal wave number kx by sin xcex8=ckx/xcfx89. A better dip response requires a longer, and hence more expensive, filter. In addition, the explicit operator method becomes more expensive and complicated in 3D compared to 2D. A third approach for affecting downward continuation with variable velocity is based on perturbation of the solution (2.3) for constant velocity. The solution (2.3) is modified according to the difference between a constant background velocity used, and the actual velocity in the subsurface. The simplest version is the split step or phase shift plus correction (PSPC) method discussed in P. Stoffa, J. Fokkema, R. Freire, and W. Kessinger. Split-step Fourier migration. Geophysics 55, 410-421, 1990, the contents of which are hereby incorporated by reference. In the PSPC method the propagation from one depth level z to the next level z+dz is obtained by applying the first term the solution (2.3) (only up going waves) with a constant background velocity c0 to             p              ~                  ~          ~                      ⁡          (                        k          x                ,                  k          y                ,        z        ,        ω            )        .
This is followed by application of a spatial Fourier transform to the x,y domain, and a multiplication of the result by the factor exp(ixcfx89dz(1/c(x,y,z)xe2x88x921/c0). The accuracy of this approach depends upon the magnitude of the difference between the values of the background velocity and the actual subsurface velocity. More accurate methods use a higher order perturbation to the constant velocity solution. These methods have been termed screen and pseudo screen propagators, and are discussed in M. De Hoop, J. Rousseau, and R. S. Wu. Generalization of the phase screen approximation for scattering of acoustic waves. Wave Motion, 31, 43-70, 2000, the contents of which are hereby incorporated by reference. The main drawback of these methods is that they are based on an approximate solution to the wave equation that becomes inaccurate in the presence of strong velocity contrasts. A partial remedy to the above-described problem is to use multiple background velocities instead of one, and at each subsurface location to select the solution corresponding to the background velocity closest in value to the actual velocity at that location.
The above solution is only partial, and methods based on the direct wave equation approach remain highly complex to calculate, although they give the most accurate results.
There is thus a widely recognized need for, and it would be highly advantageous to have, a seismic imaging method based on wave propagation but devoid of the above limitations.
In accordance with a first aspect of the present invention, a seismic imaging method of the wave equation migration type is provided by solving equation (2.2) above using an approach based on a rational Krylov expansion.
According to one aspect of the present invention there is provided a method of imaging using wave propagation data, comprising
carrying out wave migration to model wave propagation layerwise, the migration comprising:
using a krylov space expansion of an exponent of a square root operator to approximate a predetermined wave equation, and
using the approximation to model propagation of a wave between a first depth level and a second depth level, and
using the migration, together with the received data to map positions of features through which the wave propagation has occurred.
The method preferably comprises obtaining the krylov space expansion using at least one shift parameter for providing solutions of the predetermined equations in depth.
Preferably, the seismic migration is carried out on a plurality of signal traces and further comprising carrying out stacking between the traces.
In an embodiment, the stacking is carried out after the migration, in which case the data comprises M data and P data and the migration is carried out on a slice of the M data and a slice of the P data to provide M and P propagated data slices.
The method preferably comprises combining the M and P propagated data slices.
In an alternative embodiment, stacking is carried out before migration. In this case the data comprises P data and migration is carried out on a slice of the P data.
The method may utilize any one of a group of data gathering methods comprising but not exclusive to common shot data gathering, common receiver data gathering, common offset data gathering, simultaneous shot geophone data gathering, zero offset imaging, and vertical seismic profiling.
Preferably, the predetermined wave equation is the acoustic wave equation.
Preferably, the imaging is seismic imaging and the features are features of sub-surface geological structure.
Additionally or alternatively, the predetermined wave equation is the elastic vector wave equation.
Additionally or alternatively, the predetermined wave equation is a visco-elastic wave equation.
Additionally or alternatively, the approximating is carried out in time for at least one of forward migration and reverse time migration.
In one preferred embodiment, the wave propagation comprises carrying out forward simulation of wave propagation in the geological structure, and the migration comprises reverse time migration.
According to a second aspect of the present invention there is provided a method of seismic imaging comprising:
obtaining wave propagation data from a seismic source at a seismic receiver,
determining from the wave propagation data false locations of seismic features indicated in the data,
carrying out seismic migration over at least two depth levels to provide corrections to the false locations, wherein the seismic migration comprises using a krylov space expansion of an exponent of a square root operator to approximate a wave equation and therefrom to model wave propagation between the at least two depth levels,
using the corrections with the false locations to form a seismic image indicating the features.
Preferably, the wave equation is any one of a group comprising the acoustic wave equation, the elastic vector wave equation and the visco-elastic equations.
The method preferably comprises obtaining the Krylov space expansion over a depth using at least one shift parameter.
Preferably, the seismic migration is carried out on a plurality of signal traces and further comprises carrying out stacking between the traces.
In one embodiment, the stacking is carried out after the migration, in which case the data comprises M data and P data and migration is carried out on a slice of the M data and a slice of the P data to provide M and P propagated data slices.
The method preferably comprises combining the M and P propagated data slices.
Additionally or alternatively, stacking is carried out before migration. In such a case, the data comprises P data and migration is carried out on a slice of the P data.
The method non-exclusively includes utilizing any one of a group of data gathering methods comprising common shot data gathering, common receiver data gathering, common offset data gathering, simultaneous shot geophone data gathering, zero offset imaging, and vertical seismic profiling.
Preferably, the seismic migration comprises wave equation migration.
According to a third aspect of the present invention there is provided seismic imaging apparatus for receiving data propagated through sub-surface structure from a seismic source to at least one seismic receiver and forming an image therefrom, the apparatus comprising:
an initial wave analyzer for forming an initial analysis from the data including determining initial locations of seismic features indicated in the data,
a migrator, associated with the imager, for carrying out seismic migration over at least two depth levels to provide corrections to the initial locations, the migrator comprising a Krylov space expander for carrying out a Krylov space expansion of an exponent of a square root operator to approximate a predetermined wave equation, and therefrom to model a wave propagation between respective depth levels of the sub-surface structure, and
a corrector, for using the corrections with the initial locations to form a corrected mapping indicating the features, the mapping being usable for forming the seismic image.
Preferably, the predetermined wave equation is any one of a group comprising the acoustic wave equation, the elastic vector wave equation and the visco-elastic wave equations.
Preferably, the Krylov space expander is configured to obtain the Krylov space expansion over a depth using at least one shift parameter.
Preferably, the migrator is configured to carry out the seismic migration on a plurality of signal traces and the apparatus further comprises a stacker for carrying out stacking between the traces.
In an embodiment, the stacker is located after the migrator, in which case the data comprises P data and the migrator is configured to carry out the migration on a slice of the P data.
In an alternative embodiment, the stacker is located before the migrator, in which case the data comprises M data and P data and the migrator is configured to carry out the migration on a slice of the M data, and a slice of the P data to provide M and P propagated data slices.
Preferably the migrator further comprises a combiner for combining the M and P propagated data slices.
The apparatus is preferably configured to utilize at least one of a group of data gathering methods comprising common shot data gathering, common receiver data gathering, common offset data gathering, simultaneous shot geophone data gathering, zero offset imaging, and vertical seismic profiling.
Preferably, the mapping is three-dimensional mapping suitable for preparation of a three-dimensional image. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The materials, methods, and examples provided herein are illustrative only and not intended to be limiting.
Implementation of the method and system of the present invention involves performing or completing selected tasks or steps manually, automatically, or a combination thereof. Moreover, according to actual instrumentation and equipment of preferred embodiments of the method and system of the present invention, several selected steps could be implemented by hardware or by software on any operating system of any firmware or a combination thereof. For example, as hardware, selected steps of the invention could be implemented as a chip or a circuit. As software, selected steps of the invention could be implemented as a plurality of software instructions being executed by a computer using any suitable operating system. In any case, selected steps of the method and system of the invention could be described as being performed by a data processor, such as a computing platform for executing a plurality of instructions.